Predator- Prey Quasi-cycles from a Path Integral Formalism 
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The existence of beyond mean field quasi-cycle oscillations in a simple spatial model of predator 
prey interactions is derived from a path integral formalism. The results agree substantially with 
those obtained from analysis of similar models using system size expansions of the master equation. 
In all of these analyses, the discrete nature of predator prey populations and finite size effects lead 
to persistent oscillations in time, but spatial patterns fail to form. The path integral formalism 
goes beyond mean field theory and provides a focus on individual realizations of the stochastic time 
evolution of population not captured in the standard master equation approach. 
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When constructing models of biological phenomenon, 
observations of stable, periodic behavior have generally 
been taken to imply that the model will contain a sta- 
ble limit cycle. In the context of ecological modeling, 
both simple heuristic arguments and field observations 
support predator- prey oscillations in ecosystems. How- 
ever, the simple differential equation (mean field) models 
of predator prey dynamics do not exhibit limit cycles 
[l], Q . Several authors have addressed this difficulty by 
developing spatial individual level models (ILMs) that 
incorporate the stochastic effects of individualpredator- 
prey interactions as in, for example, [III These 
models yield limit cycles 6| or stochastically induced cy- 
cles dependent on space 3,0, HI- However, recent work 
on a dimensional model has shown that intrinsic noise 
without space is sufficient to generate temporal oscilla- 
tions in predator-prey populations Q • Generalization of 
this work to space shows oscillations in time, but fails to 
exhibit oscillations in space 0]. 

The purpose of the present work is to develop a mod- 
ified version of the spatial ILM of predator-prey interac- 
tions in Q and analyze the oscillatory fluctuations using 
path integral techniques. Our model includes the motion 
of both predator and prey, does not have a hard con- 
straint on the number of organisms that can be present in 
a patch and will be found to have oscillations at the global 
scale consistent with previous results Q. We map the 
master equation to a bosonic field theory [U [T3, EH E3] 
to obtain a simple derivation of coupled Langevin equa- 
tions for the fluctuations of predator-prey populations. 



DEFINITION OF THE MODEL AND MASTER 
EQUATION 



Consider a single, well-mixed patch of volume V. 
Species A is a predator for species B. We then have 



the following reactions: 



B^BB 
B^tD 
AB ^ A 
AB P2 M V AA 
A^% 



(1) 



We give the rates of the two body reactions an inverse 
V dependence, which is interpreted as the volume scaling 
of the probability in a volume V that the two organisms 
will be close enough to interact. 

The above model contains a serious defect: in the ab- 
sense of predation, the prey population diverges to infin- 
ity (in mean field). Even with predators present, this de- 
fect manifests itself through the presence of non-generic, 
initial condition dependent oscillations. To overcome this 
defect, there exist a variety of options to induce a finite 
"carrying capacity" for prey. Each option has advan- 
tages depending on the predator-prey system being de- 
scribed, though many of the predictions end up being 
generic [f3j]. One option is to restrict the total patch 
population to some number TV, including empty space 
(i.e. N A + N_r+ N E = N). This is the "urn model" 
description [lj]. In spatial models, N is often chosen 
to be I, which is equivalent to a coarse graining scheme 
which takes a patch to be the space required for one or- 
ganism. When N > 1 models are generalized to space, 
a patch is a locally well mixed area. Space is added as 
diffusion between such patches. In our model, we adopt 
the perspective that a patch is a well mixed region with 
many organisms, but do not constrain the population to 
a given N, choosing instead to obtain a finite carrying 
capacity by allowing the death rate to increase with con- 
centration. Equivalently, we could have simply included 
an intraspecies competition reaction. An advantage of 
the current approach is that it avoids nonlinear diffusive 
cross terms in spatial urn models that do not seem to 
change the dynamics substantially from versions with- 
out the cross terms Q. Additionally, urn models lead 
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to complications in the interpretation of model parame- 
ters at the mean field level and in the master equation 
due to the fact that reaction rates in urn models must be 
combined with the joint probability for drawing the reac- 
tants from the urn prior to use in the master equation or 
mean field description leading to complex combinations 
of parameters [lj]. With the soft constraint applied here, 
the reaction rates have similar, predictable meanings at 
every level of description from master equation to mean 
field. 

Formally, we include the concentration dependence of 
the death rate by noting that ua = Na/V is small 



d x (n A ) = di(0) + cn A + 0(n 2 A ), c = d'(0) > (2) 
We can now write a master equation for the patch 

d t P(m, n) = d\ (— nP(m, n) + (n + l)P(m, n + 1)) 
+c(~n 2 P(m, n) + (n + l) 2 P(m, n + 1)) 
+bi(—nP(m, n) + (n — l)P(m, n — 1)) 
+pi(-mnf(m, n) + (n + l)mP(m, n + 1)) + 
P2(—mnP(m 1 n) + (m — l)(n + l)P(m — 1, n + 1)) + 

d 2 {-mP(m, n) + (m+ l)P(m + 1, n)) (3) 

Where to denotes the number of predators, and n de- 
notes the number of prey. This master equation defines 
the time evolution of the probability distribution of pop- 
ulation states. 



MAPPING TO PATH INTEGRAL 
FORMULATION 

To analyze the predator prey dynamics, we map Eq. 
[3] to a field theory. This is done using the standard Doi 
formalism to obtain a second quantized Hamiltonian [q[ 
and bosonic coherent states to map the resulting theory 
to a path integral. For our approach and helpful reviews, 
see 

state vector 



lil li| . The mapping is achieved by introducing the 



\ip) = P(m, n)\m, n) 

m,n 

and the operator pairs a, a, b, b such that 

a\m, n) = m\m — 1, n) 
a\m, n) = \m + 1, n) 
[a, a] = 1 
b\m, n) — n\m, n — 1) 
b\m, n) = |m, n + 1) 
6, Si =1 



(4) 



Finally, all other commutators are zero. We can then 
rewrite the dynamics given by the master equation (Eq. 
[3]) as a Schrodinger like equation. 



dt\il>) = -H(a,a,b,bM 



(6) 



We now can no w sp ecify the Hamiltonian (more accu- 
rately Liouvillian [12j |) operator by multiplying the mas- 
ter equation by the state vector |m,n), summing over m 
and n, and applying the algebra of Eq. [6] to replace m 
and n by various combinations of the operators a, a and 
6, b. From this algebra, working out the structure of 
the Hamiltonian is direct and simple. As an example, we 
work out the term corresponding to prey birth explicitly 



b\ y^(-ftP(m, n) + (n — l)P(m, n — n) 
= bi J2(-bbP(m, n) + {n~- 1)P( m, n — l))|m, n) 
= —bibb\\jj) + nP(m, n)|m, n + 1) 

= -ftiS&lV) + hWb\i/>) (7) 

Other terms are treated analogously. With normal or- 
dering, this leads to the Hamiltonian 



H = bi(bb - b 2 b) + dx{bb - b) + y{b 2 b 2 - bb 2 ) 

+— i(aaS& — aab) + ^-(aabb — a 2 ab) 

+d2(aa — a) (8) 

Expectation values of functions of the random vari- 
ables to and n are given by 



(f) = (0,0\e a+b f(a,a,b,b)e 



-H(a,a,b,b)t 



|V(o)> 



(9) 



(5) 



Using bosonic coherent states, we write Eq. Has a 
path integral resulting in a Lagrangian description of the 
dynamics with generalization to space 1^, U | . Since we 
are interested in persistent oscillations around the only 
stable fixed point in the system, our choice of initial con- 
ditions is irrelevant and can be ignored. To link patches 
together for a spatial description, we define a lattice of 
patches and demand that each organism carry out a ran- 
dom walk on the lattice with given hopping probabilities 
for predator and prey. The continuum limit of a random 
walk is well known to be diffusion. We thus define diffu- 
sion rates D\ and D 2 for predator and prey respectively 
and add diffusion operators to the Lagrangian. Careful 
manipulation of the field operators leads to the same re- 
sults, provided the hopping probability for a species r 
scales as t ~ 1/a 2 where a is the lattice constant taken 
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to in the continuum limit. Then D = lim 
resulting Lagrangian density is given by 



.0 Cb T. 



The 



C = a*d t a + b*d t b - D ia *V 2 a - D 2 b*V 2 b 

+H{b,a,b,a) (10) 

With fields derived from boson operators, the La- 
grangian form of the master equation is not simply in- 
terpreted. This is because the field variables in the La- 
grangian are not simply related to the physical variables 
of population number. This proves to be the source of 
difficulties in deriving correlation functions that are phys- 
ically meaningful. To address this difficulty, we use a 
standard semi canonical Cole-Hopf transformation 17 1 
to transform the field variables to density variables 



variables. The fields z and p have a mean field value 
of due to conservation of probability [lj| • This means 
that within the Gaussian approximation, the leading or- 
der term in those fields is a small correction of order 
l/VV as above. 

To derive the mean field theory and the fluctuations, 
we then insert the rhs forms of the fields in Eq. [14] into 
the Lagrangian Eq. [13] and retain only leading and next 
to leading order, resulting in an effective Lagrangian of 
the form 



c = Wc 1 + c 2 + o(i/W) 



(15) 



Deriving each of these terms is straightforward. For 
purposes of illustration, we will carry out the expansion 
for the prey birth term explicitly 



a = ze , a = e 
b = pe-'\ I = ef 



(11) 
(12) 



This formulation has the advantage that z and p can be 
directly interpreted as the density variables for predator 
and prey respectively, while p and z generate noise terms 
at quadratic order. The transformed Lagrangian takes 
the form 



A2 



&i(v0 + vvo( 



p_ 

V 2V 1 

= b 1 (-Wp<f>~^-p V ) 



(16) 



Carrying this out for each term in the Lagrangian and 
collecting terms yields at order y/V 



C = zd t z + pd t p - D x zV 2 z - D lZ {Vz) 2 
-D 2 p(Vp) 2 - D 2P V 2 p - Ml ~ e ") 
+d 1 p(l - < ■■) - — /r( I — ( 



+ ^P(1 



r*) + ^zp(l - 



3 ) 



+d 2 z(l - e~ z ) 



(13) 



In this form, the Lagrangian has diffusive noise, and 
difficult to handle exponential terms. In the following 
section, we exploit the small parameter l/V to resolve 
these difficulties and analyse the theory. 



DERIVATION OF MEAN FIELD THEORY AND 
QUASI-CYCLES FROM LARGE V EXPANSION 

From the Lagrangian in Eq. 1131 we can proceed di- 
rectly by rewriting the fields as 



z = Vlp + VVti, p = V(f>+ \/V£ (14) 

and inserting them into the Lagrangian. These forms are 
intended to capture Gaussian fluctuations in the spirit 
of the traditional system size expansion of the master 
equation 18| while directly manipulating the population 



Ci = pd t (p + zd t if - D 1 z\7 2 tp - D 2 pV 2 (j) 
-h4>p + ditpp + cp<j) 2 + piptpcf) + P2P4"P 
—p 2 z4>ip + d 2 z<p 



(17) 



Minimizing this term provides the mean field theory. 
For V — ► oo, this minimum is exact. The Euler-Lagrange 
equations are: 



SCi 
Sz 
SCi 



= dtf — DiV 2 (p — p 2 (f>ip + d 2 tp = 
= d t <j> - D 2 V 2 (j) - b 1( j) + dx4> + c4 2 

+Plip<f> + p 2 (f)Lp = 



(18) 



These are the standard Lotka-Volterra equations gen- 
eralized to include space. They do not satisfy the crite- 
ria for pattern formation in predator-prey equations (re- 
viewed in 19]), which generically require more complex 
predation interactions. The long time dynamics relax to 
spatially uniform predator-prey populations with magni- 
tudes given by the fixed points of the ordinary differential 
equations obtained by dropping the diffusion operator in 
Eqs. [TBI above. 

At next to leading order, we fourier transform and 
switch to matrix notation, defining 



y = 



(19) 
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By simply collecting terms as in Eq. [14] we can write 
down C 2 as 



C 2 = iwy x + y Ax -y By 



The matrices are given by 



D x k 2 



-P2<f 



(Pi+P2)<p D 2 k 2 +cq 



and 



(20) 



(21) 



The power spectrum contains a nontrivial peak in uj 
corresponding to the expected temporal oscillations. The 
peak in k is at wavenumber as can be seen from the 
strictly increasing functions of k present in the spectrum. 
This rules out spatial pattern formation. These results 
are in qualitative agreement with results from expansion 
of the master equation Urn models @, H|. Additional 
work will investigate the scaling of population fluctua- 
tions near extinction transitions and in disordered envi- 
ronments. These applications are of clear ecological inter- 
est and are difficult to study with system size expansions. 
However, they can be studied using well known methods 
from field theory in the functional integral formalism. 



B = 



2(d 2 + D 1 k 2 

~P2<P<t> 



Hp 



-p 2 ip<j) 
2(fe 1 +D 2 fc 2 ) 



(22) 



We now note that the vector y is a response field 
in the Martin Siggia Rose res pon se function formal- 
ism for Langevin equations [20|, |2l|. Thus the fluctua- 
tions around mean field in the path integral are coupled 
Langevin equations. The resulting Langevin equations 
with the appropriate noise and correlations are 



— iu)x = Ax + j(uj) 



(23) 



These equations are of the same form as the equations 
reported in 0, and are easily solved using simple lin- 
ear algebra manipulations [13] 



x=-(A + iw) 1 7(w) = T>(lu) 1 j(uj) 
^x x =r l = -det(D)- 1 (D ll7l - D 1272 ) 

x 2 =e = -det(D)- 1 (D 2l7l -D 2272 ) (24) 

To obtain information from these solutions, we calcu- 
late the average power spectrum which captures oscil- 
lations but is free of phase cancellations [3]. The aver- 
age power spectrum is obtained by taking the amplitude 
squared and averaging. For predator fluctuations this 
gives 



a k + /3 k w 2 



{uj 2 



& k Y 



I> 2 



(25) 



with 



a k = Bn(k)A 2 22 + B 22 {k)A{ 2 
A = Bn(fc) 

, = D 1 k 2 {D 2 k 2 + aft) + p 2 {pi + p 2 )cj>ip > 

F = -An - A 22 
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